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20.  Abstract  (cont’d). 

A  one -dimensional  diffusion  wave  flow  routing  model,  modified  for  tailwaters,  simulates  the  important  physical  pro* 
cesses  affecting  the  flow  and  is  straightforward  to  apply.  The  model  is  based  upon  a  numerical  solution  of  the  kine¬ 
matic  wave  equation.  Ther  “modified  equation^ Hi rt,  and  von  Neumann  analyses  are  used  to  gain  insight  into  the 
stability  and  dissipative  and  dispersive  behavior  of  the  numerical  solution,  and  results  of  these  analyses  are  compared. 

A  set  of  linear  routings  is  used  to  demonstrate  the  dissipative  and  dispersive  behavior  predicted  by  the  analyses  and  to 
verify  the  accuracy  of  an  expression  that  quantifies  the  numerical  diffusion  of  the  model.  The  analyses  provide  a  basis 
for  selection  of  numerical  parameters  for  model  applications.  The  capability  and  accuracy  of  the  model  are  enhanced 
when  physical  wave  diffusion  is  balanced  by  numerical  diffusion  in  the  model.  Maintaining  the  diffusion  balance  re¬ 
quires  that  the  time  derivative  weighting  parameter  j?  be  variable  and  in  some  instances  negative.  Though  some  amount 
of  phase  error  is  introduced,  negative^ values  have  no  adverse  effect  upon  model  stability.  Field  studies  were  con¬ 
ducted  to  demonstrate  the  benefits  of  careful  model  development  and  analysis,  and  to  verify  the  diffusion  wave  model 
for  rapidly  varying  tailwater  flow.  The  bed  slope  and  roughness  characteristics  of  the  field  study  reaches  (below  Apa- 
lachia  and  Norris  Dams)  differ  greatly,  spanning  those  of  a  large  number  of  rivers  of  practical  interest.  The  accurate 
simulation  of  flow  in  both  of  these  tailwaters  attests  to  the  soundness  of  both  the  physical  basis  of  the  model  and  the 
numerical  solution  technique —The  field  studies  confirm,  for  the  extreme  case  of  rapidly  varying  flow  in  a  mildly  sloped 
river,  that  inertia  has  a  negligi^meffect  upon  unsteady  flow  waves  at  low  Froude  numbers.  Additionally,  these  studies 
verify  that  diffusion  of  short-penod  waves  in  rivers  is  generally  significant. 
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ANALYSIS  OF  A  DIFFUSION  WAVE  FLOW 
ROUTING  MODEL  WITH  APPLICATION 
TO  FLOW  IN  TAIL  WATERS 


M.G.  Ferrick,  J.  Bilmes  and  S.E.  Long 


INTRODUCTION 

Current  concerns  regarding  energy  resources  have  sparked  renewed  interest  in  hydroelectric 
power  generation.  Hydropower  is  especially  valuable  as  it  is  well  suited  to  meeting  the  peak  power 
demands  of  a  utility.  Peak  power  generation  with  hydropower  yields  flow  regimes  in  tailwater 
streams  that  are  characteristically  high  or  low  with  sharp  transitions  between  these  states.  Large 
flow  and  stage  changes  can  occur  in  a  tailwater  in  a  period  of  several  minutes.  The  duration  of  a 
zero  flow  release  can  vary  widely  with  power  demand  and  water  availability. 

Lengthy  periods  of  zero  flow  affect  the  ability  of  a  tailwater  to  maintain  a  healthy  aquatic  eco¬ 
system.  Water  temperature  and  quality  in  tailwaters  are  modified  from  those  occurring  naturally 
in  the  stream,  and  sharp  stage  transitions  can  disrupt  a  stable  river  ice  cover  in  northern  rivers. 
Accurate  knowledge  of  the  flow  regime  is  an  important  component  of  an  assessment  of  the  poten¬ 
tial  effects  of  peak  power  generation  upon  alternative  uses  of  the  stream.  In  addition,  an  under¬ 
standing  of  downstream-propagating  sharp-fronted,  large-amplitude  flow  waves  of  relatively  short 
period  is  important  because  of  their  similarity  to  dam  break  waves. 

Numerical  models  can  be  used  to  investigate  the  flow  regime  of  a  tailwater.  The  development 
of  such  a  numerical  model  has  two  basic  parts  that  we  will  address  in  this  report.  The  first  is  the 
construction  of  the  mathematical  statement  of  the  physical  processes  of  primary  interest.  The 
second  basic  part  is  the  development  and  analysis  of  the  numerical  solution  technique.  It  is  this 
second  component  of  model  development  that  often  does  not  receive  adequate  attention.  As  a 
result,  the  behavior  of  a  model  is  not  well  understood,  and  guidance  is  unavailable  concerning  pa¬ 
rameter  selection  to  achieve  optimal  accuracy  for  a  given  application  and  interpretation  of  model 
output. 

We  performed  an  analysis  of  the  dynamic  open  channel  flow  equations  to  obtain  insight  regard¬ 
ing  the  physical  processes  of  importance  in  tailwater  flow.  The  relative  magnitudes  of  wave  con¬ 
vection,  diffusion  and  dispersion  in  channels  are  represented  in  terms  of  variables  characteristic  of 
the  channel  and  the  flow.  These  processes  arc  the  result  of  contributions  of  terms  in  the  momentum 
equation,  and  when  expressed  in  nondimensional  form,  their  relative  magnitudes  indicate  appro¬ 
priate  simplifications  of  this  equation.  For  example,  the  magnitude  of  the  dimensionless  physical 
diffusion  coefficient  of  open  channel  flow  waves  is  helpful  when  considering  the  justification  for 
basing  an  analysis  upon  the  diffusion-free  kinematic  wave  theory.  In  general,  relatively  short-per¬ 
iod  waves  in  rivers  are  significantly  affected  by  diffusion.  The  analysis  also  indicates  that  inertia 
has  a  small  effect  upon  flow  waves  in  natural  channels  at  relatively  small  Fronde  numbers.  This 


conclusion  is  supported  by  the  analyses  of  Ponce  et  al.  (1978)  and  Henderson  (1%3)  for  How  in 
tailwaters,  but  contradicts  the  general  belief  that  inertia  is  important  in  rapidly  varying  flows. 

The  analyses  provided  physical  insights  that  guided  our  model  selection.  The  inertia-free  diffu¬ 
sion  wave  flow  routing  model  of  Koussis  ( 1976)  was  chosen,  and  modifications  were  made  as  neces¬ 
sary  for  application  to  tailwater  flow.  The  continuity  equation  that  forms  the  basis  of  the  model 
is  a  quasi-linear  hyperbolic  equation  for  discharge  as  a  function  of  position  and  time,  hxact  solu¬ 
tions  of  this  equation  do  not  exhibit  the  diffusion  necessary  to  simulate  wave  movement  in  most 
natural  rivers.  Through  analysis  of  the  numerical  solution,  however,  it  is  possible  to  quantify  numer¬ 
ical  diffusion  and  dispersion.  The  capability  of  the  model  is  enhanced  when  the  numerical  diffusion 
and  dispersion  error  terms  are  equated  with  the  terms  for  physical  diffusion  and  dispersion  of  (low 
waves  developed  from  the  dynamic  equations. 

The  analysis  of  a  numerical  model,  a  basic  step  in  model  development,  is  frequently  limited  to 
the  development  of  criteria  that  ensure  a  stable  solution.  Numerical  models  must  be  stable  if  the 
solution  obtained  is  to  be  meaningful.  Numerical  stability  requires  that  errors  intro-hiced  in  the 
solution  do  not  increase  in  magnitude  as  the  computation  progresses.  The  cond'.t :  quired  lor 

stability  of  a  numerical  scheme  are  frequently  known,  and  numerical  instability  f  erally  appar¬ 
ent.  The  von  Neumann  and  Hirt  analyses  have  been  used  to  develop  stability  con.  "  for  many 
numerical  models  (Roache  Id 76). 

Numerical  solutions  of  the  unsteady  open  channel  flow  equations,  however,  t;  '  '  exhibit 
errors  in  both  amplitude  and  phase  that  may  not  be  apparent  without  further  arte  >■  .umerical 
dissipation  or  diffusion  causes  the  Fourier  components  of  the  solution  and  the  erroo  io  be  damped. 
Numerical  dispersion  results  when  the  wave  celerity  of  each  wavelength  component  differs.  These 
differences  in  celerity  tend  to  modify  the  wave  form  as  the  computation  proceeds.  The  effects  of 
numerical  dissipation  and  dispersion  upon  the  accuracy  of  the  solution  are  subtle  and  difficult  to 
interpret.  An  improved  understanding  of  the  dissipative  and  dispersive  behavioi  of  the  numerical 
model  enables  the  analyst  to  minimize,  negleci  or  exploit  these  effects  to  enhance  model  accuracy 
and  to  better  interpret  computed  results.  In  our  case,  numerical  diffusion  must  be  quantified  and 
controlled  to  permit  a  model  without  physical  diffusion  to  properly  simulate  llow  that  may  he  sig¬ 
nificantly  affected  by  diffusion. 

Though  only  strictly  applicable  to  linear  equations,  the  '‘modified  equation"  (Warming  and 
Hyatt  1974)  and  von  Neumann  analyses  arc  used  to  quantify  the  dissipative  and  dispersive  behavior 
of  the  model  in  terms  of  parameters  ot  the  numerical  solution.  The  analyses  ate  complementary, 
each  having  particular  strengths.  The  Hirt  analysis  is  also  used,  and  its  model  behavior  predictions 
are  compared  with  those  of  the  other  methods.  A  set  of  linear  routings  is  used  to  demonstrate  the 
predicted  behavior  of  the  model  and  to  verify  the  adequacy  of  the  expression  for  numerical  diffu¬ 
sion  developed  in  the  modified  equation  analysis. 

Finally,  we  demonstrate  the  benefits  of  careful  model  development  and  analysts  by  comparing 
diffusion  wave  model  simulations  with  extensive  field  data  from  the  Apalachia  and  Norris  Dam 
tailwaters.  These  tailwaters  have  very  different  bed  slope  and  roughness  characteristics,  spanning 
those  of  a  large  number  of  rivers  of  practical  interest.  The  accuracy  achieved  with  the  model  in 
these  field  applications  verifies  its  generality  for  this  problem  class,  and  reinforces  the  utility  of  both 
the  nondimensiona!  statement  of  the  dynamic  equations  and  the  linear  analyses  of  numerical  solu¬ 
tion  behavior. 

PHYSICAL  DIFFUSION  AND  DISPERSION  IN  OPEN  CHANNEL  FLOW 

The  development  of  a  mathematical  statement  describing  the  important  physical  processes  of  a 
problem  relies  upon  a  clear  physical  understanding.  In  this  section  we  will  develop  a  framework  for 
obtaining  physical  insights  from  the  one-dimensional  dynamic  equations  of  flow  in  open  channels. 


Row  in  unstratified  or  weakly  stratified  reservoirs  and  in  rivers  having  a  significant  base  now  is 
generally  modeled  using  the  dynamic  equation*  S,.uiidard  numerical  solutions  of  these  equations 
fail,  however,  it  the  tlow  depth  lppiuaches  zero.  This  condition  is  common  in  tailwaters  of  dams 
used  to  generate  peak  power,  motivating  the  search  for  an  alternate  mathematical  statement. 

The  dynamic  equations  ot  tlow  in  open  channels  (St.  Venant  equations)  are  the  most  complete 
statements  of  the  laws  of  conservation  of  mass  and  momentum  is  common  use  when  the  longitu¬ 
dinal  direction  is  the  important  spatial  dimension.  The  dynamic  equations  for  a  free-flowing  river 
with  a  wide  prismatic  rectangular  channel  and  no  local  inflow  are 
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where  r  =  time 

x  =  distance  along  the  length  of  the  channel 
.1’  =  tlow  depth 
Q  =  discharge 
B  =  channel  width 
g  =  acceleration  due  to  gravity 
S0  =  slope  of  the  channel  bottom 
C=  Chezy  conveyance  coefficient. 

If  the  coefficients  of  eq  1  and  2  are  assumed  constant  at  appropriate  reference  values,  the  equa¬ 
tions  can  be  combined  and  expressed  in  terms  of  a  single  ^  'pendent  variable  yielding 
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a  hyperbolic  equation.  Equation  3  can  be  manipulated  further  to  eliminate  the  second-order  mixed 
and  temporal  derivatives  giving 

Ip  +  c  S  =  D  +  E  +  (Higher  Order  Terms) 
dr  dx  dx2  dxj 
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where  c  is  wave  celerity,  A  the  channel  cross-sectional  area  and  F  the  Froude  number,  lls/gy  • 
Kinematic  and  zero-inertia  equations  can  be  developed  as  simplified  versions  of  the  dynamic  equa¬ 
tions.  It  second  and  higher  order  terms  in  eq  4  are  neglected,  the  approximation  is  termed  kine¬ 
matic  and  is  free  of  physically  based  diffusion.  When  third  and  higher  order  terms  are  neglected, 
eq  4  contains  "ositive  physical  diffusion  and  is  j  form  of  the  zero-inertia  model.  Tracing  through 


the  development  of  eq  4  reveals  that  the  source  of  the  diffusion  term  is  primarily  the  water  surface 
slope  term  of  the  momentum  equation  (eq  2).  The  dependence  of  D  upon  the  Froude  number  re¬ 
sults  from  including  the  inertia  terms  in  the  development.  The  physical  dispersion  coefficient  E 
given  in  eq  4  is  also  positive.  Its  magnitude  varies  linearly  with  the  magnitude  of  the  diffusion  co¬ 
efficient,  and  quadratically  with  the  Froude  number.  The  existence  of  the  dispersion  term  and  the 
higher  order  terms  follow  from  the  inertia  terms. 

Equation  4  combines  the  opposing  tendencies  of  wave  diffusion  and  steepening  due  to  nonlinear 
convection.  Whitham  (1974)  showed  that  discontinuities  in  the  flow  (shock  waves)  are  not  pos¬ 
sible  if  positive  diffusion  is  present.  A  physical  justification  is  needed  for  the  application  of  kine¬ 
matic  wave  theory,  in  which  diffusion  is  neglected,  or  for  a  more  complete  description  of  the  flow. 
With  reference  discharge  Q0  and  spatial  and  temporal  increments  A x  and  At,  eq  4  is  rewritten  in 
dimensionless  form  in  terms  of  Q*  =  Q/Q0,x *  =  xf  Ax,  and  t*  =  t/At  as 
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where  Ct  is  the  Courant  number,  D*  a  dimensionless  diffusion  coefficient,  and  E*  a  dimensionless 
dispersion  coefficient.  The  magnitude  of  Ax  is  a  characteristic  of  the  wavelength.  The  value  selected 
should  provide  adequate  resolution  of  the  features  of  all  flow  waves  of  interest.  Magnitudes  of  At 
and  Q0  are  not  as  critical  as  they  appear  in  the  same  way  in  each  term  of  eq  5. 

The  magnitude  of  D*  relative  to  1  is  a  measure  of  the  importance  of  diffusion  relative  to  con¬ 
vection.  When  this  quantity  is  significantly  less  than  1 ,  convection  is  dominant  over  diffusion. 
Smooth,  steep  channels  are  therefore  good  candidates  for  shock  formation  and  successful  applica¬ 
tion  of  kinematic  wave  theory.  Other  analyses  addressing  the  use  of  the  kinematic  wave  approxi¬ 
mation  are  given  in  Menendez  and  Norscini  (1982)  and  Ponce  et  al.  (1978).  The  magnitude  of 
E*  relative  to  1  measures  the  importance  of  dispersion  relative  to  diffusion.  As  E*  is  proportional 
to  the  square  of  the  Froude  number,  its  magnitude  is  generally  much  less  than  1 .  Consequently, 
wave  dispersion  in  rivers  is  usually  negligible. 

In  natural  rivers,  flow  generally  occurs  at  small  Froude  numbers.  When  the  Froude  number  is 
significantly  less  than  1,  eq  5  simplifies  to  a  nondimensional  convective  diffusion  equation.  Cunge 
(1969)  developed  an  equivalent  dimensional  equation  by  initially  neglecting  inertia.  In  the  form 
of  eq  5,  the  dynamic  equations  indicate  that  river  flow  is  independent  of  inertia  in  many  instances. 


MODELING  APPROACH 

The  analysis  of  the  previous  section  indicates  that  an  inertia-free  simplification  of  the  dynamic 
equations  frequently  contains  all  important  physical  processes.  Flood  routing  in  rivers  and  several 
other  unsteady  open  channel  flow  problems  have  been  adequately  treated  using  simplified  zero- 
inertia  or  kinematic  wave  models.  Still,  Cunge  et  al.  (1980)  state  that  routing  methods  based  upon 
simplified  equations  may  not  be  applicable  in  situations  like  flow  in  tailwaters  where  rapid  stage 
and  discharge  variations  occur. 
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The  question  of  the  validity  of  using  simplified  routing  methods  for  rapidly  varying  tailwater 
flow  can  be  investigated  within  the  framework  provided  by  the  analyses  of  Ponce  and  Simons 
(1977),  Ponce  et  al.  (1978)  and  Henderson  (1963).  The  linear  analysis  of  Ponce  et  al.  provides 
insight  into  attenuation  and  pr  "'agation  characteristics  of  the  simplified  models  relative  to  those 
of  the  full  dynamic  model.  For  a  range  of  channel  and  flow  parameters  characteristic  of  tailwaters. 
attenuation  and  propagation  errors  resulting  from  neglecting  inertia  appear  small. 

Henderson  (1963)  compared  terms  of  the  momentum  equation  in  an  order-of-magnitude  anal¬ 
ysis  for  a  wide  rectangular  channel.  The  acceleration  terms  were  of  the  same  order  of  magnitude 
and  were  related  to  the  water  surface  slope  term  as 


acceleration  term 
by/bx 


0(F2). 


As  the  Froude  number  of  the  (low  in  natural  rivers  is  typically  much  less  than  I .  the  acceleration 
terms  are  small  relative  to  the  water  surface  slope  term.  Henderson  also  compared  the  water  sur¬ 
face  and  bottom  slope  terms,  obtaining 

bvjbx  a  bqjbi 

So  “  qV*S 05/3 

where  q  is  the  discharge  per  unit  width.  For  flood  flows  in  natural  streams  the  magnitude  of  the 
water  surface  slope  term  is  often  small  relative  to  the  bottom  slope  term.  In  tailwaters,  where 
sudden  large -magnitude  flow  changes  occur,  bq/br  can  be  large,  and  the  water  surface  slope  term 
can  be  expected  to  make  an  important  contribution  to  the  momentum  balance.  The  importance 
of  the  water  surface  slope  term  is  further  enhanced  if  the  river  being  modeled  has  a  small  bottom 
slope. 

Tire  analyses  support  the  use  of  an  inertia-free  model  for  tailwater  flow  routing.  Our  analysis 
and  that  of  Ponce  et  al.  (1978)  indicate  the  importance  of  diffusion,  and  the  Henderson  analysis 
yields  the  related  result  that  water  surface  slope  is  an  important  part  of  the  momentum  balance. 
The  zero-inertia  model  satisfies  these  requirements.  The  local  and  convective  acceleration  terms 
of  the  full  momentum  equation  are  neglected  in  this  model.  The  abbreviated  set  of  equations  in¬ 
cludes  the  water  surface  slope  term  which  permits  wave  attenuation  or  diffusion,  and  can  be  solved 
without  difficulty  as  the  flow  depth  becomes  small.  Commonly  in  tailwaters.  however,  the  water 
surface  slope  changes  abruptly  with  the  passage  of  the  toe  of  a  wave  front,  and  a  number  of  fronts 
may  exist  simultaneously.  As  a  result,  the  complications  of  accurately  locating  each  front  and 
evaluating  the  water  surface  slope  in  the  vicinity  of  a  front  must  be  introduced  into  the  model. 

The  kinematic  wave  approximation  is  obtained  by  neglecting  the  water  surface  slope  term  in 
addition  to  the  acceleration  terms  of  the  momentum  equation.  This  additional  simplification  does 
not  permit  wave  attenuation  and  can  be  a  serious  limitation  if  the  river  being  modeled  is  long  or 
mildly  sloped.  Numerical  solutions  of  the  kinematic  wave  equation,  however,  frequently  exhrbit 
wave  attenuation  resulting  from  the  solution  technique  (Cunge  1969,  Smith  1980).  This  numer¬ 
ical  diffusion  can  be  used  to  mimic  the  physical  diffusion  occurring  in  the  channel  and  overcome 
the  inherent  model  limitation  to  diffusion -free  flows. 

Simple  model  formulation  makes  application  of  a  kinematic  wave  based  routing  technique  to 
flow  in  tailwater  streams  attractive.  Models  based  upon  the  kinematic  wave  approximation,  in¬ 
cluding  the  well-known  Muskingum  model,  determine  flow  at  a  point  independent  of  downstream 
influences.  These  models  do  not  require  a  downstream  boundary  condition  and  cannot  account 
for  the  influence  of  downstream  controls  upon  the  flow.  The  ability  of  this  group  of  models  to 
correctly  represent  flow  in  rivers  with  long  pooled  reaches  is  therefore  suspect.  Smith  (1980) 
found  that  a  variable  weighting  factor  in  the  numerical  scheme,  corresponding  to  a  variable  dif¬ 
fusion  coefficient,  allowed  successful  application  of  kinematic  wave  based  models  to  flood  routing 


through  flat,  ponded  reaches.  He  did  not,  however,  resolve  the  importance  of  backwater  effects 
upon  rapidly  varying  flows  found  in  tailwaters. 

In  conclusion,  the  tailwater  flow  model  formulation  should  include  variable  wave  attenuation 
and  the  contribution  of  the  water  surface  slope  to  the  momentum  balance.  As  a  physically  mean¬ 
ingful  downstream  boundary  is  not  generally  available  in  tailwaters,  models  not  requiring  a  down¬ 
stream  boundary  condition  are  most  readily  applied.  The  diffusion  wave  model  that  is  described, 
analyzed  and  tested  in  this  report  satisfies  these  requirements.  The  model  is  an  adaptation  for 
tailwaters  (Ferrick  1980)  of  the  Hood  routing  model  of  Koussis  (1976).  Weinmann  and  Laurenson 
(1979)  related  this  formulation  to  other  simplified  routing  methods. 


DESCRIPTION  OF  THE  DIFFUSION  WAVE  FLOW  ROUTING  MODEL 


The  diffusion  wave  model  for  flow  in  tailwaters  is  a  kinematic  wave  based  approach  that  differs 
from  most  simplified  routing  methods  in  that  both  stage  and  discharge  are  computed  at  each  point 
in  the  numerical  grid.  The  conservation  of  mass  equation  is  solved  numerically  for  discharge.  Ex¬ 
act  solutions  of  the  continuity  equation  deform  but  do  not  damp  with  time,  so  that  diffusion, 
which  is  necessary  for  representation  of  important  physical  processes  in  tailwater  Hows,  is  not  present 
explicitly.  The  numerical  solution  of  this  equation,  however,  is  adjusted  in  the  model  to  require 
numerical  diffusion  to  mimic  physical  diffusion.  An  equation  for  the  stream  rating  curve,  which 
includes  the  zero-inertia  form  of  the  momentum  equation,  models  river  siage.  The  inclusion  of  the 
water  surface  slope  term  in  this  relationship  generates  a  looped  rating  curve  and  provides  an  im¬ 
proved  estimate  of  wave  celerity.  The  pair  of  equations  arc  coupled  through  the  wave  celerity  and 
must  be  solved  simultaneously. 

The  conservation  of  mass  equation  for  flow  in  open  channels  can  he  w  i  it  ten  as 


d Q 

dx 
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where 


_  dQ  _  dx 
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Q  =  discharge, 

q(  =  local  inflow  per  unit  length  of  the  channel 
A  =  cross-sectional  area  of  the  channel 
c  =  wave  celerity 

Equation  7,  which  forms  the  basis  of  the  model,  is  a  first-order  hyperbolic  equation.  This  type  of 
equation  is  advantageous  for  modeling  tailwater  flow  because  a  downstream  boundary  condition 
need  not  be  specified.  For  this  same  reason,  however,  backwater  effects  upon  the  flow  cannot  be 
taken  into  account. 

The  “method  of  lines”  approach  is  used  to  obtain  a  solution  of  eq  7.  In  this  method  the  spatial 
derivative  is  approximated  with  a  finite  difference  expression,  but  the  dependent  variable  remains 
continuous  in  the  time  domain.  The  partial  differential  equation  for  conservation  of  mass  is  there¬ 
by  reduced  to  an  ordinary  differential  equation.  With  the  approximations 
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in  which  j  is  an  index  corresponding  to  the  spatial  location  x  =  /Ax,  and  8  is  a  parameter  of  the 
numerical  solution,  eq  7  is  rewritten  as 

6J+1  (/)  +  *  Gj+1  «  =  a  Q^+aAx  q{t)  -  j  Qt(t)  (10) 


where 


<c> 

a=  Ax(l  -0)  ' 

Q  is  the  derivative  of  discharge  with  respect  to  time,  and  <c>  is  an  averaged  celerity  in  space  and 
time  over  the  cell  of  the  computational  mesh  where  the  equation  is  applied.  Equation  10  is  a  linear 
first-order  ordinary  differential  equation.  Assuming,  for  a  small  time  increment,  that  the  variation 
of  Gj(f)  linear,  that  the  values  of  the  coefficients  can  be  estimated,  and  that  the  local  inflow  is 
constant,  the  solution  of  eq  10  is 


CjTl'  =  ci  +  C4 

where 

Ct  =  1  -  a 
C2=a-0 

C3  =  0-  exp  (£j) 

C4  =<7  Ax  (1-0) 
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The  local  Courant  number  of  the  numerical  grid  cell  Cr  expresses  the  ratio  of  physical  to  numer¬ 
ical  wave  celerity,  and  m  is  the  temporal  index  of  the  computational  grid,  t  =  mAt. 

The  variable  physical  diffusion  of  the  flow  is  represented  by  tuning  the  numerical  diffusion  in¬ 
herent  in  the  modfel.  An  expression  for  the  model  parameter  0  that  results  from  enforcing  the 
physical/numerical  diffusion  balance  is 


X  = 


<£>. 


BS0  <c>  Ax 


(1  -F2/4) 


02) 


where  discharge  and  wave  celerity  are  averaged  over  the  numerical  grid  cell,  and  B  is  the  width  of 
the  channel.  A  necessary  condition  for  stability  of  the  model  is  that  either  0  <  0.5  or  0  <  1.0 
and  Cr  >  1.0.  Equation  12  and  the  stability  criteria  for  the  model  will  be  developed  subsequently. 
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For  small  values  of  the  Froude  and  Courant  numbers,  eq  12  is  equivalent  to  the  expression  devel¬ 
oped  by  Cunge  (1969)  for  the  weighting  factor  in  the  Muskingum  model: 


.  <Q>  N  . 

BS0<c>  Ax/ 


(13) 


Neglecting  the  inertia  of  the  flow  and  the  momentum  contribution  of  the  local  inflow,  the 
zero-inertia  conservation  of  momentum  equation  for  prismatic  channels  is 
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in  which  Sf  is  the  slope  of  the  energy  grade  line. 

The  zero-inertia  momentum  equation  inserted  into  Manning's  equation  yields  an  expression  for 
the  stream  rating  curve, 


where  R  is  the  hydraulic  radius,  n  the  Mannings  roughness  coefficient,  and  C’m  a  constant  that 
is  dependent  upon  the  system  of  units.  To  obtain  a  relationship  that  is  consistent  with  the  down¬ 
stream-progressing  discharge  calculation  (eq  1 1 ),  the  spatial  derivative  in  eq  I  5  is  replaced  with  a 
quantity  determined  at  a  point.  If  the  energy'  slope  is  adequately  large,  arguments  from  kinematic 
wave  theory  can  be  used,  and  eq  15  can  be  rewritten  in  a  form  of  the  “Jones  formula"  as 
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With  a  finite  difference  approximation  of  the  time  derivative  in  eq  16  and  the  assumption  of  a 
wide  rectangular  channel,  the  flow  depth  is  determined  as 
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The  accuracy  of  the  flow  depth  calculation  using  eq  1 7  depends  upon  the  validity  of  the  kine¬ 
matic  wave  approximation.  In  response  to  decreasing  flow  in  a  river  reach  having  a  small  bottom 
slope,  the  denominator  of  eq  17  may  decrease  more  quickly  than  the  numerator,  causing  the  cal¬ 
culated  stage  to  increase.  This  unphysical  result  signals  the  need  lor  an  alternate  equation  for  the 
calculation  of  flow  depth.  If  the  definition  of  celerity  is  approximated  as 

dA  '  8  Ay 

then  an  alternate  equation  for  How  depth  can  be  found: 

y=y0  +  (<? '&<>)■  <18) 

The  remaining  unknown  to  be  determined  is  wave  celerity.  As  discussed  by  Henderson  (1963), 
the  celerity  of  a  kinematic  or  zero-inertia  flood  wave  is  related  to  the  flow  velocity  by  a  multiplier 
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that  depends  upon  the  channel  shape  and  energy  slope  model  used.  The  wave  celerity  of  steep- 
fronted  tailwater  releases  is  also  dependent  upon  the  flow  depth  on  either  side  of  the  wave  front. 
In  the  extreme  case  of  a  rapid  flow  release  to  a  previously  dry  channel,  the  celerity  of  the  front 
must  equal  the  velocity  of  the  flow  immediately  behind  it.  The  monoclinal  rising  wave  is  a  trans¬ 
lator  wave  of  stable  form.  If  rising  tailwater  releases  are  presumed  to  be  monoclinal  waves  during 
passage  through  a  reach  Ax,  then  an  expression  for  wave  celerity  can  be  obtained.  For  a  wide  rec¬ 
tangular  channel,  with  the  Chezy  equation  used  to  describe  the  energy  slope,  the  expression  for 
wave  celerity  is 


p-cw^l 


(19) 


For  slowly  rising  hydrographs,  eq  19  yields  the  familiar  result  for  celerity  of  a  flood  wave,  c  =  1 .5 
V,  On  the  recession  side  of  the  hydrograph,  the  wave  celerity  used  in  the  model  is  the  value  for 
flood  waves. 

The  water  surface  slope  is  largely  responsible  for  diffusion  of  flow  waves  and  the  existence  of 
a  looped  rating  curve  in  rivers.  In  the  model,  however,  diffusion  is  generated  numerically  with  the 
6  parameter  and  does  not  depend  directly  upon  the  water  surface  slope  term.  The  depth  corre¬ 
sponding  to  a  given  discharge  varies  depending  upon  the  evaluation  of  the  water  surface  slope  (eq 
1 7),  and  is  lower  during  the  rising  limb  of  a  hydrograph  than  during  the  falling  portion  of  the  hy¬ 
drograph.  This  in  turn  affects  the  modeled  discharge  and  wave  celerity  through  eq  1 1  and  19,  re¬ 
spectively. 

To  advance  the  routing  one  time  step,  eq  6  is  solved  at  each  grid  point  using  values  of  Q  and  c 
computed  at  the  previous  time  step  to  evaluate  the  coefficients.  The  values  of  flow  depth  and 
celerity  are  then  updated  with  eq  17  through  19.  The  updated  values  are  used  to  reevaluate  the 
coefficients  of  eq  1 1 ,  and  the  iteration  repeats.  A  converged  solution  is  reached  when  each  of  the 
flow  values  computed  at  successive  iterations  agrees  within  a  set  tolerance.  In  simulations  per¬ 
formed  at  small  Courant  numbers  a  good  initial  estimate  of  the  solution  is  available  from  the  pre¬ 
vious  time  step,  and  convergence  is  generally  rapid. 

The  parameter  0  varies  locally  with  the  flow  and  is  continually  updated  during  the  iterative 
process.  Values  of  the  wave  celerity  and  the  corresponding  Courant  number  vary  along  the  chan¬ 
nel.  The  model  time  step  is  determined  from  a  limitation  imposed  upon  the  maximum  Courant 
number.  The  Courant  number  limitation,  required  to  preserve  accuracy,  has  a  numerical  value 
that  is  dependent  upon  flow  conditions. 


ANALYSIS  OF  THE  NUMERICAL  MODEL 

The  development  of  a  numerical  model  should  include  a  thorough  analysis  of  the  solution  tech¬ 
nique.  We  will  analyze  the  diffusion  wave  model  with  several  linear  techniques  in  this  section  of 
the  report.  The  difference  expression  for  the  continuity  equation  used  in  our  model  was  com¬ 
pared  by  Koussis  (1980)  with  the  Muskingum  model,  a  more  widely  used  method  for  routing  flow 
in  open  channels.  The  difference  representation  of  this  equation  is  repeated  here  for  the  case  of 
no  local  inflow, 

*  0  -o)CP  +  l  +(a-0)fij"  +  /3G£,  (20) 

0  =  e*p(iTe) 
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where  m  and  j  are  indices  of  the  temporal  and  spatial  location  on  the  computational  grid,  respec¬ 
tively. 

The  diffusion  wave  model  as  applied  by  Weinmann  and  Laurenson  ( 1 979)  and  as  described 
above  is  nonlinear  with  coefficients  that  are  dependent  upon  the  local  flow  conditions.  The  cal¬ 
culation  for  Q  is  explicit,  progressing  in  the  downstream  direction  from  a  known  flow  rate  at  the 
upstream  boundary  and  given  initial  conditions.  At  each  new  time  level  the  model  iterates  until 
the  change  in  computed  flow  at  each  grid  point  is  below  a  specified  tolerance. 

Modified  equation  and  Hirt  analyses  of  diffusion  wave  model 

Stability,  damping  and  dispersion  characteristics  of  a  difference  approximation  to  a  partial  dif¬ 
ferential  equation  can  be  investigated  with  the  modified  equation  analysis  of  Warming  and  Hyett 
and  with  the  Hirt  analysis.  The  two  analyses  follow  the  same  basic  steps  with  one  important  ex¬ 
ception. 

Neglecting  roundoff  error,  the  modified  equation  represents  the  actual  partial  differential  equa¬ 
tion  solved  when  a  numerical  solution  is  obtained  from  a  difference  equation.  To  obtain  the  mod¬ 
ified  equation  for  the  difference  scheme  of  eq  20,  each  term  is  expanded  in  a  Taylor  series  about 
Q™ .  Upon  simplification,  the  resulting  equation  is 

dQf  dQ  |  cAx  3 ]Q_  +  At  3^+  cAt  d2Q  +  cAx2  <PQ  +  Ax^  d3Q 
to  C  Zx*  2  dx2  2  dr2  <l~fl  dxdr  6  dx3  2 a  3*23, 

,  Ax  At  b3Q  ,  At2  d3Q  ,  L  „  ,  T 

+  — —  - —  +  — —  — -  +  (Higher  Order  Terms)  =0  (21 ) 

2a  dxdt2  6  df3 

The  modified  equation  has  an  infinite  number  of  terms.  Terms  which  appear  in  the  modified  equa¬ 
tion  but  are  missing  from  the  original  differential  equation  represent  a  type  of  truncation  error. 

Properties  of  a  difference  scheme  can  be  found  by  examining  a  truncated  version  of  the  modi¬ 
fied  equation.  The  time  derivatives  higher  than  first  order  and  the  mixed  derivatives  are  eliminated 
from  eq  21  to  obtain  an  equation  that  is  amenable  to  physical  interpretation.  Kven-order  spatial 
derivatives  in  the  recast  equation  correspond  to  dissipative  effects  and  odd-order  spatial  derivatives 
reveal  dispersive  properties  of  the  model.  In  the  Hirt  analysis,  the  governing  differential  equation 
is  used  to  simplify  eq  21 .  A  solution  of  the  original  differential  equation  will  not,  in  general,  sat¬ 
isfy  the  difference  equation.  Therefore,  in  the  modified  equation  approach,  eq  21  itself  is  differ¬ 
entiated  and  used  in  the  simplifying  process.  The  coefficients  are  assumed  to  be  constant  in  both 
analyses.  Differences  between  the  two  procedures  for  the  diffusion  wave  model  analysis  appear 
in  the  coefficients  of  third  and  higher  order  spatial  derivatives. 

Following  the  modified  equation  approach,  eq  21  becomes 


+  (Higher  Order  Terms) 

to  dx  dx2  3*3 


D  = 

E  = 


-I  -2C?  (l- 


1  "0  O-0): 


Jl 


Following  Hirt’s  analysis,  the  analogous  expression  for  I:  is 


.  _  cAx2 
H“  6 


(22) 


(23) 
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If  third  and  higher  order  terms  are  ignored,  eq  22  is  a  linear  convective  diffusion  equation.  The 
Fourier  components  of  the  continuum  solution  of  this  equation  are 

fik(x-Ct)-Dk2t 

where  k  is  the  wave  number  (2rr/wave  length)  of  the  component,  i  =  \f~\,  and  Qk  is  the  amplitude 
of  the  component  of  wave  number,  k.  The  modified  equation  of  the  diffusion  wave  model  (eq  22) 
can  be  rewritten  in  the  form 


f  ■  t  <«<*'> 

p=0 


d2P+lQ 

QX1P*1 


+ 


£  M2 P) 

p=l 


&PQ 
a x*p 


The  form  of  the  solution  of  eq  24  is 


Q(x,t)^^a*‘b>t  <*kx 


where 


<*=  £  (-l)P*JPM2p),  b=  £(-1)PA:jp+1  M2p  +  1). 


p=i 


p  =  0 


(24) 


(25) 


As  waves  of  large  wave  number  cannot  be  resolved  on  the  numerical  grid,  waves  having  small 
wave  numbers  are  of  primary  importance  in  the  numerical  solution.  For  these  waves  the  exponent 
a  of  eq  25  can  be  approximated  as 


a  =«  -k2D  (26) 

where  D  is  the  diffusion  coefficient  defined  in  eq  22.  A  dimensionless  numerical  diffusion  coef¬ 
ficient  D*  =  D/(cAx )  of  the  diffusion  wave  model  is  given  in  Figure  1  as  a  function  of  the  Courant 
number  for  various  values  of  9.  A  positive  diffusion  coefficient  in  the  modified  equation  is  neces¬ 
sary  for  positive  damping  in  the  model  and  a  stable  numerical  solution.  A  stable  solution  is  ob 
tained  if  either  9  <.  0.5  or  Cr  >  1 .0  and  9  <  1.0.  Stability  does  not  impose  a  restriction  upon  the 
maximum  value  of  the  Courant  number,  which  presents  the  possibility  of  using  large  time  steps 
in  the  model.  Numerical  dissipation  increases  as  9  decreases  and  as  Cr  increases.  In  the  low  Cour¬ 
ant  number  range,  <  0.5,  damping  does  not  vary  greatly  with  the  Courant  number. 

Components  of  all  wave  numbers  of  the  continuum  solution  of  the  convective  diffusion  equa¬ 
tion  are  advected  at  c.  After  an  increment  of  time  A r  the  components  have  each  undergone  a 
phase  angle  change,  4>c,  of 

4>c  =  -dtA/  =  - Cr  (*Ax)  =  -  Cr  y  .  (27) 

For  the  difference  approximation,  celerity  is  a  function  of  wave  number  (eq  25).  The  phase  angle 
change  of  the  numerical  solution  in  time  At  is 

om 

<hN  =  bAt  -  k  fj(\)At  ■¥  At  £  (-1)P  **P+1  M2P+1).  (28) 

p=i 


The  ratio  of  the  numerical  to  continuum  phase  shifts  yields  an  expression  for  the  relative  propa¬ 
gation  speed  of  each  Fourier  component  per  time  increment 


(29) 


*,  =  V* c  - 1  - l-  E  HP  *2p  *2P  +  »  • 

p=i 


Values  of  4>r  greater  than  1  indicate  that  the  numerical  solution  component  of  wave  number  k  will 
have  a  celerity  greater  than  that  of  the  continuum  solution.  The  converse  is  true  for  values  of  'h, 
less  than  1.  Since  small  wave  numbers  are  of  primary  importance  in  the  numerical  solution,  eq  29 
can  be  approximated  as 


.=  1  + 


k 2  m(3) 


+  0(y4)  *  1  +  Y 


[3f.(rr ('-iVtfbp)]' 


(30) 


Equation  30  is  plotted  as  a  function  of  the  Courant  number  for  selected  values  of  8  in  Figures  2-4 
for  wavelengths  of  24  Ax',  12  Ax  and  6  Ax,  respectively.  For  these  wavelengths,  the  phase  angle 
of  the  numerical  solution  predominantly  lags  that  of  the  continuum  solution.  The  discrepancy  is 
largest  for  the  shorter  wavelengths,  smaller  values  of  6  and  larger  values  of  the  Courant  number. 
For  the  shorter  wavelengths  and  Courant  numbers  greater  than  1,  the  phase  angle  of  the  numerical 
solution  varies  strongly  with  the  Courant  number.  The  numerical  solution  is  most  likely  to  exhibit 
leading  phase  angles  for  waves  of  short  wavelength  with  6  approximately  0.5. 


Figure  1 .  Dimension  lew  numerical 
diffusion  as  a  function  of  Courant 
number  for  various  values  of  the  pa¬ 
rameter  8,  Values  of  D*  are  based 
upon  the  modified  equation  analysis. 
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Figure  2.  Ratio  of  numerical  to  con¬ 
tinuum  phase  shifts  in  time  At  for 
24 -Ax  wavelengths  as  a  function  of 
Courant  number  and  various  values 
of  6.  One  set  of  curves  is  based  upon 
the  modified  equation  analysis  and 
one  set  upon  the  von  Neumann  analysis. 
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Figure  3.  Ratio  of  numerical  to 
continuum  phase  shifts  in  time 
At  for  1 2- Ax  wavelengths  as  a 
function  of  Courant  number 
and  various  values  of  6.  One 
set  of  curves  is  based  upon  the 
modified  equation  analysis 
and  one  set  upon  the  von 
Neumann  analysis. 


Figure  4.  Ratio  of  numerical  to 
continuum  phase  shifts  in  time 
At  for  6-  Ax  wavelengths  as  a 
function  of  Courant  number 
and  various  values  of  0.  One 
set  of  curves  is  based  upon  the 
modified  equation  analysis 
and  one  set  upon  the  von 
Neumann  analysis. 
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Figure  3.  Ratio  of  numerical  to 
continuum  phase  shifts  in  time 
At  for  24-Ax,  1 2-Ax  and  ft-Ax 
wavelengths  as  a  function  of 
Courant  number  and  i arums 
values  of  0.  The  curves  are 
based  upon  the  llirt  analysis. 


An  analogous  equation  for  the  ratio  of  the  numerical  to  continuum  phase  shifts  based  upon  the 
Hint  analysis  can  be  written  by  substitution  of  eq  23  for  p  (3)  of  eq  30.  yielding 


3  cr  / 

—r)  -  1  1 

'  *  ->)l 

r  \ 

1  -0/  r  ' 

a  -  a  /I 

Equation  31  is  plotted  in  Figure  5  as  a  function  of  the  Courant  number  for  the  same  wavelengths 
and  0  values  used  with  eq  30.  This  series  of  curves  projects  quite  different  dispersion  behavior 
than  the  modified  equation  analysis.  At  a  Courant  number  of  1 .  the  Hirt  analysis  projects  zero 
dispersion  independent  of  wavelength  and  6.  Lagging  phase  angles  are  projected  at  higher  Courant 
numbers,  and  all  phase  angles  are  leading  at  lower  Courant  numbers.  Shorter  wavelengths  and 
lower  values  of  6  are  projected  to  have  larger  phase  errors. 


von  Neumann  analysis  of  the  diffusion  wave  model 

As  a  result  of  the  truncation  of  terms  in  eq  26,  the  modified  equation  analysis  does  not  provide 
information  regarding  the  diffusive  nature  of  the  short  wavelength  Fourier  components  of  the  so¬ 
lution  or  of  the  error.  These  components  are,  at  times,  responsible  for  instability  of  a  numerical 
solution.  Similarly,  the  approximation  in  eq  30  limits  the  application  of  the  modified-equation- 
based  phase  relationship  to  waves  of  relatively  small  wave  number.  In  addition  to  its  usual  role  of 
providing  numerical  stability  criteria,  the  von  Neumann  analysis  can  be  used  to  identify  the  relative 
damping  and  dispersion  of  Fourier  components  of  all  wave  numbers. 

The  evolution  of  the  numerical  solution  in  a  time  step  Ar  is  considered  in  this  approach.  The 
coefficients  a  and  0ofeq  20  are  assumed  to  be  constant.  The  solution  ofeq  20  can  then  be  written 
as  a  Fourier  series: 


<r  =£e?  e'n 

k 


(32) 


where  is  an  amplitude  function  at  time  in  At  of  the  Fourier  component  of  wave  number  k. 
Each  term  of  the  difference  equation  is  replaced  by  its  A-th  Fourier  component.  The  spatial  do¬ 
main  is  assumed  to  be  infinite,  and  boundary  influences  are  not  considered.  The  decay  or  ampli¬ 
fication  of  each  component  is  evaluated  to  investigate  stability  and  damping  of  the  numerical 
scheme  by  forming  the  ratio  of  the  amplitude  functions  at  two  successive  times,  nAl  and  (n  +  1  )A/. 
Performing  these  operations  upon  eq  20  yields 

Q™  a  +  (en  -  1  ) 


in  which  the  complex  number  rk  is  termed  the  amplification  factor.  A  rreeessary  and  sufficient 
condition  for  stability  of  the  solution  is  that  the  modulus  of  rk  be  less  than  or  equal  to  1  for  all 
integer  values  of  A;  (Richtmyer  1957).  An  expression  for  the  square  of  the  modulus  of  rk  which 
follows  from  eq  33  is 


a4  +  2a2  (0+  1  -a)  (0+  1 )  ( I  -cos  7)  +  20  |2ji  ( 1  -  a)  +  a  (a-  2  )|  ( I  —cos  >)2  -  2  a 2  p  sin2  7 
a2  +  4  a2  ( I  -  a)  ( 1  -cos  7)  +  4  ( 1  -  a)2  ( I  -cos  7)2 


(34) 


The  square  of  the  modulus  of  rk  must  also  remain  less  dian  or  equal  to  1  for  all  integer  values  of  k 
for  a  stable  numerical  solution.  When  |rk  P  is  equal  to  1  for  a  component  of  wave  number  k,  the 
numerical  scheme  is  termed  conservative  or  neutrally  stable.  Smaller  values  of  |rk  P  correspond  to 
larger  inherent  dissipation  of  the  numerical  scheme. 

Long  wavelengths  are  resolved  over  a  large  number  of  grid  points.  For  these  waves,  7  is  small 
and  the  following  substitutions  can  be  made  in  eq  34: 

sin2  7==72 

l  -  cos 7  *  72/2 

(1  -  COS7)2  =*  0  (35) 


and  eq  34  can  then  be  simplified  to  yield 


(36) 


I  h  a2  +  [(1  -q)  +  fl0?-a)b2 
k  a2 +  2(1 -a)?2 

The  solution  will  be  stable  for  long  period  waves  if 

0(0-<*)-(l  -<*)<0.  (37) 

This  inequality  is  satisfied  if  either  0  <  0.5  or  Ct  >  1 .0  and  8  <  1 .0.  These  criteria  are  in  agreement 
with  those  based  upon  the  modified  equation  analysis.  Figure  6  presents  |rR|2  as  a  function  of 
Courant  number  for  two  wavelengths  (24Ax,  1 2Ax)  and  selected  values  of  8  (0.9, 0.5, 0.0,  -1 .0). 
Numerical  damping  is  sensitive  to  Cr  and  0,  as  was  noted  in  the  modified  equation  analysis,  and 
exhibits  the  same  trends.  For  long  period  waves  having  equivalent  6  and  Cr,  the  shorter  wave¬ 
lengths  are  more  hi^tly  damped. 

As  mentioned  above,  the  von  Neumann  approach  permits  analysis  of  short  period  waves. 
Assuming  y  is  equal  to  jt/2,  eq  34  becomes 

L  |2  _  o2  -2a3  (1  +<3)+  2a2  (1  +  g)2  -4afl(l  +  (3)  ■+  4(32 
k  a4  -4  a3  +8  a2  -8a  +  4 


Figure  6.  Square  of  the  modulus 
of  the  amplification  factor  for 
24 -Ax  and  1 2-  Ax  wavelengths 
as  a  function  of  Courant  number 
and  various  values  of  0. 


Figure  7.  Square  of  the  modulus 
of  the  amplification  factor  for 
4-  Ax  and  2- Ax  wavelengths  as 
a  function  of  Courant  number 
and  various  values  of  8. 
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Imposing  the  stability  restriction  on  the  square  of  the  modulus  of  rk  requires  that 

aJ  -  a2  (0  +  3)  +  2  a(0  +  2)  -  2  (0  +  1 )  <  0  .  (39) 

The  stability  limits  developed  for  the  long  period  waves  also  satisfy  the  requirement  imposed  by 
eq  39  for  a  short  period  wave.  The  |rk  P  is  presented  as  a  function  of  Courant  number  in  Figure  7 
for  values  of  7  equal  to  both  n/2  and  tr,  ti  *  shortest  wavelength  resolvable  on  the  numerical  grid. 
Short  period  wave  damping  increases  with  Courant  number  and  is  a  much  stronger  function  of 
Courant  number  than  for  the  long  period  waves.  Also,  damping  generally  increases  as  wavelength 
decreases.  Wave  propagation  with  a  Courant  number  of  1 .0  and  0  of  0.9  is  undamped  for  both 
short  and  long  period  waves.  No  limitation  upon  the  minimum  value  of  6  is  necessary  for  stability 
of  the  model. 

The  dispersive  properties  of  the  numerical  scheme  can  also  be  investigated  using  the  von  Neu¬ 
mann  analysis.  For  small  angles,  the  measure  of  an  angle  in  radians  is  approximately  equal  to  the 
sine  of  the  angle.  The  phase  angle  of  the  numerical  solution  at  time  At  is  then 

•mag  (rk) 

*N  =  •  (40) 

For  the  diffusion  wave  model,  eq  40  can  be  written  as 


|‘1>N  I  =  sin  y  •  C  - - 1 

C2  sin2  y  +  I  +  2/1  +  A  2  I 


Cr(/32  +2/3Cf-  1)(1  -cos7) 

A  O-0)J 

The  ratio  of  the  numerical  to  continuum  phase  angle  is 

4,  =  JELL  f _ • _ 1  v' 

^  |_Cr2  sin2  7  +  1  +  14  + /42J 

Again,  propagation  in  the  numerical  model  matches  that  in  the  continuum  solution  of  the  convec¬ 
tive  diffusion  equation  if  this  ratio  is  1  for  a  given  wavelength.  Equation  42  is  plotted  in  Figures 
2-4  for  wavelengths  of  24Ax,  12  Ax  and  6Ax  respectively.  Though  the  values  of  4>r  are  not  gen¬ 
erally  equal  to  those  obtained  using  the  modified  equation  analysis,  all  of  the  trends  are  in  agree¬ 
ment.  Unlike  the  modified  equation  approach,  however,  no  restriction  upon  7  was  used  in  the 
development  of  the  equation  for  4^.  Therefore,  for  the  shorter  wavelengths  (Fig.  4),  the  von 
Neumann  analysis  is  likely  to  give  a  better  estimate  of  the  phase  variation  in  the  model.  Numeri¬ 
cal  dispersion  is  typically  the  largest  for  short  period  waves  or  short  period  Fourier  components  of 
a  wave.  Equation  42  reveals  that  the  shortest  wave  resolvable  on  the  numerical  grid  is  stationary. 


LINEAR  CASE  STUDIES 

A  set  of  linear  case  studies  is  presented  to  demonstrate  the  utility  of  the  modified  equation  and 
von  Neumann  analyses  in  representing  model  behavior.  Figure  8  presents  half  sine  waves,  of  wave¬ 
length  8 Ax  and  16Ax,  which  will  serve  as  initial  conditions  for  the  studies  which  follow.  In  each 
case,  the  wave  celerity  is  held  constant  at  3  ft/s  (0.91  m/s),  independent  of  the  flow.  If  in  addition 
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Figure  8.  Half  sine  waves  of  wavelength  8-  and  16- Ax 
that  serve  as  initial  conditions  for  the  linear  case  studies. 

to  constant  celerity,  a  constant  diffusion  coefficient  is  assumed,  then  the  solution  of  the  initial 
value  problem  posed  by  the  linear  convective  diffusion  equation  and  an  initial  condition  of  Figure 
8  can  be  written  as  a  Fourier  series. 

Initially,  a  key  question  will  be  addressed  regarding  the  adequacy  of  the  expression  given  by  the 
modified  equation  analysis  for  the  numerical  diffusion  inherent  in  the  model.  A  07-term  Fourier 
series  solution  of  the  convective-diffusion  equation  with  a  diffusion  coefficient  calculated  from  eq 
22  was  obtained.  This  analytical  solution  is  compared  in  Figure  9  with  the  corresponding  numer¬ 
ical  solution  for  each  wavelength  after  the  center  of  the  wave  has  been  advected  6  miles  (9.7  km) 
downstream.  Cases  presented  were  projected  in  the  analyses  to  have  minimal  dispersion  and  a 
range  of  diffusion.  The  damping  evidenced  in  the  numerical  and  analytical  solutions  are  in  agree¬ 
ment  in  all  cases.  The  cases  with  8  equal  to  0.9  and  Cr  of  1 .0  exhibited  essentially  pure  convection 
No  damping  or  phase  error  of  any  Fourier  component  of  the  solution  was  projected  for  this  case 
by  the  von  Neumann  analysis.  A  small  amount  of  dispersion,  evidenced  by  a  slight  lag  in  the  nu¬ 
merical  solution,  was  present  in  the  cases  with  0  of  0.0  and  Cr  of  0.1 .  This  slight  lagging  phase 
error  was  projected  by  both  the  von  Neumann  and  modified  equation  analyses  (Figs.  2  and  3).  By 
contrast,  the  HiTt  analysis  projected  a  leading  phase  error  for  these  same  cases  (Fig.  5). 

Several  cases  are  presented  to  demonstrate  the  projected  diffusive  and  dispersive  behavior  of  the 
model.  Figure  10  contains  numerical  and  analytical  solution  comparisons  for  cases  having  a  con¬ 
stant  Courant  number  and  selected  values  of  0,  after  the  center  of  the  wave  has  been  advected  6 
miles  (9.7  km)  downstream.  The  leading  phase  error  of  the  longer  wavelength  components  and 
the  more  extreme  leading  phase  error  of  the  short  wavelength  components  are  evident  for  the  case 
in  which  8  is  0.5.  This  behavior  was  projected  in  Figures  3  and  4.  Minimal  damping  of  the  short 
wavelength  components  (Fig.  7)  is  the  other  condition  necessary  for  development  of  leading  short 
period  waves.  A  small  amount  of  phase  lag  occurs  as  projected  for  the  case  with  9  equal  to  0.0. 

In  this  case,  damping  is  large,  especially  of  the  short  wavelength  components.  Larger  damping  and 
greater  phase  lag  occur,  as  expected,  when  0  is  set  at  -1 .0.  The  Hirt  analysis  projects  leading  phase 
errors  for  all  9  values  when  the  Courant  number  is  small.  In  addition,  leading  errors  are  projected 
to  be  more  extreme  at  small  values  of  0(<  0.0). 

For  the  cases  presented  in  Figure  1 1 , 6  is  held  constant  as  Courant  number  is  varied.  As  pro¬ 
jected  by  die  modified  equation  and  von  Neumann  analyses,  leading  phase  errors  occur  for  Cour¬ 
ant  number  0.1 ,  lagging  phase  errors  occur  for  Courant  number  4.0,  and  a  minimal  phase  error  is 
observed  for  Courant  number  1.0.  Numerical  diffusion  is  again  well  represented  by  the  expression 
developed  using  the  modified  equation  approach. 
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Figure  9.  Comparison  of  numerical  and  Fourier  series 
solutions  for  8-  and  16-Ax  wavelengths  after  the  center 
of  the  wave  has  propagated  6  miles  (9.7  km)  downstream. 
Cases  shown  were  projected  to  have  minimal  dispersion 
and  a  range  of  diffusion  by  the  modified  equation  and 
von  Neumann  analyses. 


4  8 

Distance  (miles) 


Figure  10.  Comparison  of  numerical  and  Fourier  series 
solutions  for  the  8- Ax  wavelength  and  a  fixed  value  of 
Cf  =  0. 1  after  the  center  of  the  wave  has  propagated  6 
miles  ( 9.7  km)  downstream. 


In  practice,  cost  considerations  generally  preclude  the  use  of  an  extremely  fine  computational 
grid.  An  assessment  of  the  adequacy  of  a  given  mesh  relative  to  the  flow  features  of  interest  in  the 
prototype  is  an  important  part  of  model  testing.  The  linear  analyses  have  shown  that  numerical 
damping  and  phase  errors  increase  with  spatial  grid  size.  To  investigate  these  effects,  the  short 
period  wave  of  Figure  8  was  resolved  on  a  coarse  grid  containing  one-half  the  number  of  points 
used  in  Figure  8.  Cases  having  parameters  identical  to  those  presented  in  Figure  10  were  repeated 
using  this  coarse  grid  and  are  presented  in  Figure  12.  As  diffusion  increases  with  the  spatial  mesh 
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Figure  II.  Comparison  of  numerical  and  Fourier  ser¬ 
ies  solutions  for  the  16- wavelength  and  a  fixed 
value  of  6  =  0.5  after  the  center  of  the  wave  has  prop¬ 
agated  6  miles  (9.7  km)  downstream. 
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Figure  12.  Comparison  of  numerical  and  Fourier  series 
solutions  for  the  shorter  (8-Yx)  wavelength  resolved  on 
a  coarse  grid  after  the  center  of  the  wave  has  propagated 
6  miles  (9. 7  km)  downstream.  The  wave  had  a  4-  Ax  in¬ 
itial  wavelength  relative  to  the  coarse  grid. 


dimension  Ax,  damping  is  greater  in  the  coarse  mesh  simulation  for  given  values  of  8  and  the  Cow- 
ant  number.  Large  diffusion  cases  are  well  behaved;  the  wavelength  of  the  flow  quickly  increases 
and  numerical  and  analytical  solutions  correspond  as  before.  The  lightly  damped  case  (O  =  0.1 , 

6  =  0.5)  largely  retains  its  short  wavelength  and  exhibits  the  same  tendencies  as  its  finer  meshed 
counterpart.  The  peak  flow,  however,  does  not  correspond  to  the  analytical  solution  and  the 
amplitude  of  the  leading  waves  has  increased.  The  wave  period  of  interest  is  short  enough,  rela¬ 
tive  to  the  numerical  grid,  that  a  more  accurate  estimate  of  numerical  diffusion  is  required.  By 
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retaining  more  terms  of  the  modified  equation,  a  hitter  order  estimate  of  diffusion  could  be  made. 
It  can  be  observed  in  Figure  10,  however,  that  a  refined  mesh  tends  to  correct  the  diffusion  problem 
and  yields  decreased  amplitudes  of  the  leading  waves.  Adequate  accuracy  of  the  numerical  solution 
on  a  coarse  grid  requires  the  presence  of  a  large  amount  of  diffusion.  This  requirement  could  pre¬ 
clude  or  impair  the  analysis  of  physically  important  cases. 


ACCURACY  CONSIDERATIONS  OF  THE  NUMERICAL  SOLUTION 


The  order  of  accuracy  of  the  difference  scheme  is  defined,  in  the  modified  equation  analysis,  as 
the  power  of  the  computational  mesh  dimension  in  the  coefficient  of  the  lowest  order  error  term. 
If  D  of  eq  22  is  physically  based,  the  numerical  scheme  is  accurate  to  second  order.  Physically 
based  values  of  D  and  E  yield  a  third-order  numerical  solution. 

Cunge  (1969)  and  Smith  (1980)  equated  coefficients  of  physical  and  numerical  diffusion  to  en¬ 
hance  the  capabilities  of  the  Muskingum  flow  routing  model.  Koussis  (1976)  also  followed  this 
approach  with  the  present  diffusion  wave  model.  In  equating  the  diffusion  coefficients,  an  error 
term  in  the  numerical  solution  of  the  continuity  equation  is  used  to  extend  the  physical  basis  of 
the  model.  An  equation  for  the  numerical  weighting  factor  0,  in  terms  of  Courant  number  and 
physical  parameters  of  the  flow  and  channel,  is  obtained  by  equating  the  diffusion  coefficients  of 
eq  4  and  22: 


0  =  1  + 


Cr 


In 


/ 1  +  X  -  Cr 
\1  +  A  +  Cr 


) 


(43) 


where 


The  value  of  0  is  a  function  of  position  and  time.  For  flow  at  relatively  low  Froude  numbers  com¬ 
mon  in  rivers,  eq  43  reduces  to  the  equation  given  by  Koussis.  As  evidenced  in  Figures  1, 6  and  7, 
model  damping  is  dependent  upon  the  Courant  number  and  0.  As  the  Courant  number  approaches 
1 ,  the  balance  between  physical  and  numerical  diffusion  is  maintained  by  allowing  0  to  increase  to 
compensate  for  increased  damping  as  a  result  of  the  Courant  number.  As  the  Courant  number  is 
increased  further,  0  calculated  from  eq  43  may  exceed  the  upper  bound  for  model  stability.  In 
this  instance,  small  values  of  diffusion  cannot  be  attained  in  the  model.  If  the  diffusion  balance  is 
not  enforced,  the  accuracy  of  the  model  will  be  degraded. 

Strupczewski  and  Kundzewicz  (1980)  and  Dooge  et  al.  (1982)  found  in  analyses  of  the  Musk¬ 
ingum  model  that  negative  values  of  the  weighting  parameter  are  due  to  short  model  reach  lengths. 
The  same  conclusion  can  be  drawn  from  the  value  of  0  expressed  in  eq  13  and  43.  If  negative 
values  occur,  0  can  no  longer  be  considered  a  weighting  parameter.  Instead,  it  is  a  parameter  used 
to  control  or  tune  the  diffusion  of  the  numerical  model. 

In  the  interest  of  further  enhancing  the  physical  basis  of  the  diffusion  wave  model,  equating 
the  coefficients  of  the  physical  and  numerical  dispersion  terms  of  eq  4  and  22,  respectively,  would 
be  desirable.  This  balance  would  provide  improved  model  phase  accuracy.  However,  the  model 
does  not  contain  a  free  variable  in  addition  to  the  0  parameter  with  which  to  enforce  balanced 
dispersion.  Given  0  of  about  0.2S  and  Cr  <  0.S,  the  dispersion  term  of  the  modified  equation  is 
approximately  zero.  Larger  values  of  0  correspond  to  a  positive  numerical  dispersion  coefficient, 
and  conversely,  smaller  0  values  indicate  a  negative  coefficient.  The  physical  dispersion  coefficient 
of  eq  4  is  always  positive.  To  properly  simulate  physical  dispersion  requires  that  the  shorter  wave¬ 
length  components  propagate  more  quickly  than  the  longer  wavelength  components,  if  the 
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dispersion  coefficient  in  the  model  is  negative,  just  the  opposite  will  occur.  Minimization  of  the 
imbalance  is  possible  for  many  cases  of  interest  with  judicious  selection  of  reach  length.  An  ade¬ 
quate  description  of  sharp-fronted  waves  of  relatively  short  period  that  are  encountered  in  tail- 
waters  may,  however,  dictate  fine  resolution  of  computational  mesh.  Negative  values  of  6,  which 
may  then  be  required  for  proper  diffusion  in  the  model,  will  introduce  some  amount  of  lagging 
phase  error  in  the  simulation  that  will  be  most  severe  for  the  shorter  wavelengths.  Assessing  the 
relative  magnitudes  of  the  dimensionless  coefficients  in  eq  5  reveals  that  the  magnitude  of  the  dis¬ 
persion  term  is  small  relative  to  that  of  the  diffusion  term.  Therefore,  imprecise  modeling  of  dis¬ 
persion  is  probably  not  a  serious  limitation  In  any  case,  the  insights  obtained  from  these  analyses 
remain  valuable  as  the  lagging  phase  error  can  be  anticipated,  improving  the  interpretation  of  model 
output. 

FIELD  STUDIES 

Two  extensive  field  tests  were  conducted  to  verify  the  diffusion  wave  model.  One  test  was  con¬ 
ducted  in  the  Hiwassee  River  below  Apalachia  Dam  and  the  other  in  the  Clinch  River  below  Norris 
Dam.  Apalachia  Dam  is  situated  at  Hiwassee  River  mile  (HRM)  66  near  the  southern  end  of  the 
Tennessee/North  Carolina  border.  At  the  dam,  the  river  has  a  drainage  area  of  1018  miles2  (2636 
km2).  In  the  13-mile  (20.9-km)  study  reach  below  the  dam  the  river  bed  has  a  steep  slope,  dropping 
over  350  ft  (106.7  m).  The  bed  is  extremely  rough  with  large  boulders  and  trees  in  the  channel, 
creating  roughness  elements  that  are  typically  on  the  order  of  the  flow  depth.  The  pools  located 
in  the  reach  are  relatively  short.  The  bulk  of  the  flow  in  this  reach  is  normally  diverted  from  the 
dam  through  a  conduit  to  the  Apalachia  powerhouse  located  near  the  end  of  the  reach.  Therefore 
flow  occurs  only  during  floods  and  as  a  result  of  local  inflows  from  a  drainage  area  totaling  1 18 
miles2  (306  km2). 

Norris  Dam  is  situated  at  Clinch  River  mile  (CRM)  79.7  near  Oak  Ridge,  Tennessee.  The  drain¬ 
age  area  of  the  river  at  Norris  Dam  is  2912  miles2  (7542  km2).  The  physical  characteristics  of  the 
Norris  tailwater  are  significantly  different  from  those  of  the  Apalachia  tailwater.  The  tailwater 
has  a  relatively  mild  bed  slope,  dropping  about  25  ft  (7.6  m)  in  the  1 3-mile  (20.9-km)  study  reach 
immediately  below  the  dam.  Relative  to  that  of  the  Apalachia  Dam,  a  much  greater  percentage  of 
the  Norris  tailwater  is  pooled  at  low  flow.  An  individual  pool  located  near  the  center  of  the  stud> 
reach  is  over  2.5  miles  (4.0  km)  in  length  and  has  a  bed  slope  of  only  0.00012.  In  general,  rough¬ 
ness  elements  on  the  Clinch  River  reach  are  much  smaller  than  those  on  the  Hiwassee  River  reach. 

In  the  pools,  most  of  the  roughness  elements  are  submerged  even  during  lengthy  zero  flow  release 
periods  from  Norris  Dam. 

The  physical  characteristics  of  the  Hiwassee  River  below  Apalachia  Dam  and  the  Clinch  River 
below  Norris  Dam  span  those  of  a  large  number  of  streams.  Both  rivers  have  typical  alternating 
pool-riffle  structures.  The  study  reaches  were  located  immediately  below  each  dam  where  the 
features  of  the  tailwater  hydrograph  are  the  sharpest  and  the  effect  of  diffusion  upon  the  flow  re¬ 
leases  is  most  pronounced.  As  the  diffusion  of  a  flow  wave  in  a  wide  rectangular  channel  is  in¬ 
versely  proportional  to  the  bed  slope  (eq  4),  much  greater  diffusion  of  the  flow  releases  was  ex¬ 
pected  in  the  Norris  tailwater  relative  to  that  in  the  Apalachia  tailwater. 

Apalachia  Dam  tailwater 

A  field  investigation  of  flow  in  the  Apalachia  Dam  tailwater  was  conducted  by  personnel  with 
the  Tennessee  Valley  Authority  (TV A)  on  22-23  March  1979.  The  hydrograph  during  this  study 
was  produced  with  sluice  gates  at  the  dam  and  is  given  in  Figure  13.  River  stage  was  recorded  at 
HRM  62.8,  59.0,  56.9  and  53.0  during  the  test.  The  channel  shape  was  approximated  in  the  model 
as  rectangular  throughout  the  tailwater.  Channel  width  and  slope  were  obtained  from  USGS 
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Figure  13.  Apakchia  Dam  flow  releases.  22-23  March  1 9  79. 

quadrangles.  The  tailwater  channel  width  varied  between  1 50  and  440  ft  (46  and  1 34  m),  averaging 
280  ft  (85  m),  and  the  channel  bed  slope  varied  between  0.0027  and  0.0084.  The  Manning’s 
roughness  coefficient  was  estimated  on  the  basis  of  field  observations  and  was  adjusted  during 
model  calibration.  Calibrated  roughness  coefficients  ranged  between  0.04  and  0.07,  averaging 
0.066.  Lateral  inflow  was  initially  neglected. 

In  addition  to  physical  parameters  characterizing  the  reach,  model  application  requires  a  selec¬ 
tion  of  numerical  grid  parameters.  Both  the  linear  analyses  and  past  experience  with  model  appli¬ 
cation  provide  guidance  to  initiate  the  selection  process.  We  performed  further  studies  of  the 
laboratory  test  cases  reported  in  Ferrick  (1980)  that  have  shown  that  mass  balance  and  wave  propaga¬ 
tion  speed  errors  in  a  model  simulation  are  reduced  with  improved  spatial  grid  refinement  and  by 
decreasing  the  maximum  Courant  number.  Linear  analysis  of  the  model  has  revealed  that  damping 
and  phase  errors  occur  when  the  Ax  chosen  is  large  relative  to  wavelengths  being  routed.  For  these 
reasons,  a  spatial  grid  size  of  5280  ft  (1609  m)  and  a  maximum  Courant  number  of  1 .0  were  chosen 
for  application  of  the  model  to  the  Apalachia  tailwater.  A  limitation  upon  the  minimum  value  of 
the  parameter  0  of  0.0  was  imposed  in  the  model  to  prevent  lagging  phase  errors  that  occur  with 
negative  0  values. 

A  comparison  of  measured  and  computed  stage  is  presented  in  Figure  14.  The  propagation  of 
the  1000-ft3/s  (28.3-m3/s)  release  in  the  model  lags  the  data  by  a  time  that  increases  with  distance 
downstream.  The  timing  and  magnitude  of  the  other  releases  are  accurately  represented  in  the 
model.  Due  to  the  absence  of  the  powerhouse  discharge,  the  measured  and  computed  stages  are 
not  in  agreement  at  HRM  53.  The  diffusion  of  the  waves  during  passage  through  the  study  reach 
is  not  large  (Fig.  1 5),  and  the  computed  rating  curves  at  the  four  stage  measurement  locations  were 
not  strongly  looped.  These  observations  indicate  that  the  water  surface  slope  term  retained  in  the 
momentum  balance  may  not  be  important.  The  water  surface  slope  increases  the  effective  channel 
slope  at  the  wave  front  and  causes  an  increase  in  the  speed  of  wave  propagation.  The  effect  of 
neglecting  this  term  upon  the  modeled  flow  was  to  lag  the  arrival  of  the  rising  limb  of  each  hydro¬ 
graph.  Though  the  shape  of  the  computed  hydrograph  was  not  significantly  altered,  discounting 
the  water  surface  slope  degraded  the  wave  timing  agreement  between  the  model  and  the  data. 

A  number  of  sensitivity  studies  were  conducted  in  which  the  estimated  input  parameters,  width 
and  roughness,  were  varied  in  an  effort  to  improve  the  agreement  of  the  modeled  and  measured 
propagation  of  the  1000-ft3/s  (28  J-m3/s)  wave.  Increased  channel  roughness  caused  a  reduction 
in  wave  speed,  an  increase  in  steady-state  stage  for  a  given  flow,  and  an  extended  duration  of  the 
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Figure  14.  Measured  and  computed  stage  at  several 
locations  on  the  Apalachia  tailwater.  A  maximum 
Courant  number  of  1 .0,  spatial  grid  resolution  of 
5280 ft  (1609  m),  and  variable  6  limited  to  pos¬ 
itive  values  were  conditions  of  the  numerical 
simulation. 


Figure  15.  Computed  discharge  at  several  locations  on 
the  Apalachia  tailwater.  A  maximum  Courant  number 
of  1.0,  spatial  grid  resolution  of 5280  ft  (1609)  m),  and 
variable  0  limited  to  positive  values  were  conditions  of 
the  numerical  simulation. 
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increased  stage.  Increased  channel  width  caused  a  decrease  in  the  magnitude  of  stage  changes,  a 
slowing  of  wave  movement  and  a  decreased  period  of  increased  flow.  The  overall  agreement  of 
model  and  prototype,  however,  was  not  improved. 

A  physically  justifiable  development  that  produced  the  correct  propagation  speed  of  the  1000- 
ft3/s  wave  without  significantly  affecting  the  larger  waves  was  the  inclusion  of  local  inflows.  Local 
inflows,  measured  at  20  ft3/s-mile  (0.35  m3/s-km)  for  the  3.5-mile  (5.6-km)  reach  nearest  the  dam, 
were  assumed  to  be  representative  for  the  tailwater.  In  addition,  the  known  discharge  at  the  power¬ 
house  was  included  as  a  local  inflow.  Figure  14  presents  a  revised  comparison  between  the  model 
and  prototype.  At  HRM  63  the  only  discrepancy  concerns  the  stage  increase  during  passage  of  the 
1000-ft3/s  (28.3-m3/s)  flow.  The  effective  channel  width  at  low  flow  is  less  than  at  higher  flow, 
but  due  to  the  assumption  of  a  rectangular  channel,  measured  stage  at  low  flow  is  higher  than  the 
computed  stage.  At  HRM  59  the  averaged  channel  parameters  of  the  model  accurately  describe 
the  local  channel  characteristics  except  at  times  when  only  local  inflows  are  present.  Sensitivity 
studies  indicated  that  the  remaining  discrepancies  in  the  arrival  times  of  the  1000-ft3/s  (28.3-m3/s) 
wave  at  HRM  59  and  57  can  be  attributed  to  a  lack  of  detailed  inflow  data.  The  river  is  pooled  at 
low  flow  at  mile  57,  and  as  a  result,  the  computed  river  stage  at  low  flow  at  this  location  falls  be¬ 
low  the  measured  stage.  Timing  and  magnitude  of  the  larger  flow  arrivals  are  again  accurately  rep¬ 
resented  in  the  model.  Further  downstream,  at  HRM  53,  model  and  prototype  data  agree  closely, 
both  in  timing  and  magnitude  for  all  the  releases. 

The  flow  releases  in  the  Apalachia  tailwater  are  lightly  damped  (Fig.  14-15),  and  experience 
with  the  model  has  revealed  that  this  condition  increases  the  tendency  to  compute  an  unphysical 
discharge  value  immediately  preceding  a  sudden  increase  or  decrease  in  discharge.  This  error  was 
suppressed  in  the  model  by  retaining  the  discharge  computed  at  the  previous  time  step.  Numerical 
experiments  were  conducted  to  verify  that  this  mechanism  did  not  adversely  affect  the  solution. 
Further  experiments  have  shown  that  the  tendency  to  compute  an  unphysical  discharge  is  reduced 
by  decreasing  the  maximum  Courant  number,  or  equivalently,  the  temporal  grid  resolution. 

Dynamic  waves  of  measurable  size  traveling  ahead  of  the  main  flow  and  upstream  wave  move¬ 
ment  resulting  from  wave  reflections  were  not  observed  during  the  test.  Together  with  these  ob¬ 
servations,  the  accurate  simulation  of  wave  propagation  and  damping  of  the  rapidly  varying  flows 
generated  for  the  Apalachia  tailwater  study  confirm  the  negligible  influence  of  the  inertia  terms 
for  flow  in  steeply  sloping  tailwaters.  The  ability  of  the  model  to  simulate  flow  in  mildly  sloping 
tailwaters,  where  the  inertia  terms  are  larger  relative  to  other  terms  of  the  momentum  balance, 
however,  is  more  important  practically  and  a  more  severe  test  case. 

From  the  hydrographs  presented  in  Figure  1 5  it  is  noted  that  the  model  accurately  conserved 
mass.  The  limitation  imposed  upon  the  6  parameter,  6  >  0.0,  permitted  the  numerical  diffusion 
in  the  model  to  maintain  the  balance  with  the  physical  diffusion  occurring  in  the  stream  and  did 
not  impair  the  accuracy  of  the  simulation.  The  Manning’s  roughness  model  for  the  energy  slope 
was  adequate  to  characterize  the  effect  of  the  large-scale  roughness  elements  upon  the  flow.  The 
physical  justification  for  omission  of  backwater  effects  in  the  model  cannot  yet  be  judged.  It  is 
possible  that  the  pooled  reaches  in  the  Apalachia  tailwater  were  short  enough  to  have  a  negligible 
effect  upon  the  flow  or  that  the  addition  of  local  inflow  compensated  for  the  effect. 

Norris  Dam  tailwater 

A  160-hour  controlled  release  test  (Fig.  16)  was  performed  in  the  Norris  Dam  tailwater  by  per¬ 
sonnel  with  the  TVA  on  1-7  July  1980,  during  which  the  variation  of  river  stage  was  continuously 
recorded  at  CRM  78.85,  76.1 , 73.6,  71 .4,  and  67.3.  The  long  pool  in  the  study  reach  was  isolated 
by  the  placement  of  two  of  the  recording  gages  at  CRM  76.1  and  73.6.  The  channel  shape  was 
assumed  in  the  model  to  be  rectangular  throughout  the  study  reach.  Prior  to  the  test  TVA  per¬ 
sonnel  surveyed  the  tailwater  for  channel  width  and  bed  slope.  Measured  widths  ranged  between 
260  and  540  ft  (79  and  165  m),  averaging  380  ft  (1 1 6  m),  and  bed  slopes  varied  between  0.0001 2 
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Figure  1 6.  Norris  Dam  flow  releases,  1-7  July  1 980. 
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Figure  1 7.  Hydrographs  at  downstream  extent  of  Norris  tailwater  study  reach  computed  with  constant  0 
0.5,  spatial  grid  resolution  of  2640  ft  (805  m),  and  maximum  Courant  numbers  of  0.25  and  1.0. 


and  0.00130,  Local  inflows  were  small  during  the  test,  averaging  about  2  ft3/s-mile  (0.035  m3/s- 
km),  with  point  inflows  of  19  ft3  (0.54  m3/s)  from  Coal  Creek,  the  largest  tributary,  and  40  ft3/s 
(1.13  m3/s)  of  leakage  past  Norris  Dam.  Calibrated  roughness  coefficients  that  were  used  in  the 
model  simulations  ranged  between  0.01 5  and  0.035,  averaging  0.026. 

The  Norris  tailwater  model  required  a  relatively  fine  spatial  resolution  of  2640  ft  (805  m)  to 
provide  stage  and  discharge  information  that  was  adequate  to  assess  the  effect  of  alternative  flow 
release  strategies  upon  the  tailwater  fishery,  and  to  provide  resolution  necessary  for  accurate  rout¬ 
ing  of  short  period  releases.  The  linear  analysis  revelaed  that  the  dissipative  and  dispersive  char¬ 
acteristics  of  the  model  are  sensitive  to  the  Courant  number  and  the  6  parameter,  in  addition  to 
the  spatial  grid  resolution.  To  achieve  optimal  model  accuracy,  the  limitations  to  be  imposed  upon 
the  maximum  value  of  the  Courant  number  and  the  minimum  value  of  the  parameter  0  were  sys¬ 
tematically  addressed.  Constant  6  simulations  for  various  values  of  the  maximum  Courant  number 
and  equal  Courant  number  simulations  with  different  minimum  0  limitations  were  run  to  supple¬ 
ment  the  information  supplied  in  the  linear  analysis. 

Figure  17  presents  computed  flow  at  CRM  67.2  with  a  constant  0  value  of  0.5  and  maximum 
Courant  numbers  of  0.25  and  1.0,  respectively.  The  simulation  at  the  higher  Courant  number  ex¬ 
hibited  greater  model  damping  and  lagging  of  the  waves,  relative  to  the  lower  Courant  number 
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Figure  18.  Hydrographs  at  downstream  extent  of  Minis  tailwater  study  reach  computed  with  constant  0  =0.0,  spatial 
grid  resolution  of  2640  ft  (805  m),  and  maximum  Courant  numbers  of  0.25  and  1.0. 


simulation.  This  effect  was  especially  pronounced  for  short  period  waves.  Model  damping  in  the 
small  Courant  number  simulation  was  minimal.  Based  upon  the  linear  model  analysis,  all  of  these 
tendencies  were  expected.  Both  simulations  created  mass,  and  the  smaller  Courant  number  case 
has  a  mass  conservation  error  of  28%  at  CRM  67.2. 

A  comparison  of  the  simulations  with  a  constant  0  value  of  0.0,  given  in  Figure  1 8,  shows  that 
wave  propagation  was  slightly  faster  and  wave  attenuation  was  apparently  greater  for  the  simula¬ 
tion  with  a  maximum  Courant  number  of  0.25  than  for  the  Courant  number  1 .0  simulation.  The 
contradiction  with  the  linear  theory  concerning  wave  attenuation  is  actually  an  effect  of  improved 
mass  conservation  which  is  seen  as  a  quicker  arrival  of  the  tail  of  the  wave.  The  bulk  of  the  dif¬ 
ference  between  the  two  simulations  occurred  in  the  long,  nearly  flat,  pooled  reach  between  CRM 

76.2  and  73.7.  Again,  both  simulations  created  mass,  but  the  error  was  only  4%  for  the  smaller 
Courant  number  case.  Measured  and  computed  stages  with  a  constant  6  value  of  0.0  and  a  maxi¬ 
mum  Courant  number  of  0.25  agree  reasonably  well  at  all  gages,  although  modeled  damping  is 
generally  less  than  that  in  the  prototype.  This  result  indicates  that  a  variable  weighting  factor  is 
not  essential  for  reasonably  accurate  modeling  of  flat,  ponded  river  reaches. 

Comparing  the  simulations  of  like  Courant  number  in  Figures  1 7  and  18  reveals  model  tenden¬ 
cies  due  to  varying  the  6  parameter  that  were  predicted  in  the  linear  analysis.  The  effect  upon 
model  damping  of  varying  0  is  much  greater  than  the  effect  of  varying  the  Courant  number.  For 
the  cases  in  which  9  is  set  at  0.0,  the  slopes  of  wave  fronts  are  less  steep,  corresponding  to  increased 
diffusion  relative  to  the  cases  where  9  is  0.5.  Also,  the  arrivals  of  the  short  wavelengths  at  CRM 

67.2  are  lagged  in  the  6  equal  0.0  cases,  relative  to  the  9  equal  0.5  cases. 

Further  reduction  of  the  Courant  number  did  not  affect  the  computed  hydrograph  for  any  of 
the  constant  6  runs  at  the  downstream  end  of  the  study  reach.  Larger  Courant  number  (>  2.0) 
simulations  were  attempted,  but  the  iteration  did  not  converge  to  a  solution  at  the  initial  abrupt 
flow  increase  of  the  inflow  hydrograph.  This  convergence  problem  was  not  resolved  because  of 
generally  poor  model  accuracy  at  large  Courant  numbers. 

The  linear  analysis  of  the  numerical  model  indicated  that  negative  values  of  9,  required  to 
maintain  the  physical-numerical  diffusion  balance  when  fine  numerical  grids  are  used,  introduce 
a  small  lagging  phase  error  in  the  modeled  results.  Computed  hydrographs  for  simulations  having 
either  a  0m)n  of  0.0  or  no  limitation  upon  0min  are  presented  in  Figure  19  for  CRM  73.7  and  67.2. 
Figure  19  reveals  that  the  timing  of  wave  arrival  is  not  greatly  affected  by  the  0mtn  limitation, 
but  wave  damping  and  mass  balance  are  extremely  sensitive  to  the  limitation.  Limiting  0  to 
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Figure  l  9.  Hydrographs  at  downstream  extent  of  the  long  pool  and  downstream  extent  of  Norris  tailwater 
study  reach  computed  with  maximum  Cowant  number  of  0.25,  spatial  grid  resolution  of  2640  ft  (805  m), 
and  variable  0  either  limited  to  positive  values  or  allowing  negative  values  ■ 


positive  values  results  in  an  increase  in  mass  of  14%  at  the  downstream  end  of  the  study  reach. 
Relaxing  the  0  limitation  to  be  greater  than  -1 .0  reduces  the  increase  of  mass  to  8%.  Removing 
the  limitation  altogether  yields  a  4%  decrease  of  mass  at  the  downstream  end  of  the  reach. 

Only  minimal  differences  between  the  simulations  existed  upstream  of  the  2.5-mile  (4.0-km) 
pooled  reach  where  bed  slope  is  relatively  large  and  the  0min  limitation  is  imposed  infrequently. 
Downstream  of  the  pool  at  CRM  73.7,  the  flow  and  stage  differences  between  the  simulations  are 
pronounced  and  continue  to  increase  to  the  downstream  extent  of  the  study  reach.  Much  larger 
wave  diffusion  is  expected  and  is  evident  in  the  results  of  the  simulation  without  a  minimum  8 
limitation.  Limiting  8  to  positive  values  did  not  produce  attenuation  of  the  modeled  stages  that 
was  adequate  to  reproduce  the  prototype  stage  measurements. 

A  comparison  of  measured  and  computed  stages  without  a  0min  limitation  at  five  locations  on 
the  tailwater  is  given  in  Figure  20.  The  stage  measurement  locations  do  not  coincide  exactly  with 
the  modeled  sections  and  proper  interpretation  of  the  offsets  between  measured  and  computed 
stage  requires  that  these  differences  be  considered.  Overall  agreement  between  the  model  and  the 
data  on  wave  timing,  damping  and  shape  is  excellent.  At  the  upstream-most  gage  (CRM  78.85) 
the  model  (CRM  78.7)  reproduces  the  shape  and  timing  of  all  the  releases.  An  error  in  the  chart 
speed  of  the  stage  recorder  beginning  at  hour  100  leads  to  an  increasing  timing  discrepancy  toward 
the  end  of  the  test. 

hi  Figure  20  the  gage  at  the  head  of  the  long  pool,  CRM  76.1 ,  is  compared  to  the  model  output 
at  CRM  76.2.  The  timing  of  the  wave  arrivals  is  well  represented  in  the  model.  The  modeled  stages 
for  the  smaller  waves  are  slightly  less  damped  than  in  the  prototype.  In  the  pool,  almost  all  rough¬ 
ness  elements  were  submerged  at  low  flow.  Small,  but  physically  reasonable,  values  of  Manning's 
roughness,  0.01 5,  were  required  to  attain  the  proper  depths  of  flow  in  the  model  for  the  pooled 


27 


»  •*» 


Stage 


T ime  (hrs) 

Figure  20.  Measured  and  computed  stage  at  several  locations  on  the  Norris  tailwater. 
A  maximum  Courant  number  of  0.25,  spatial  grid  resolution  of  2640  ft  (805  m ),  and 
variable  d  without  a  lower  bound  were  conditions  of  the  numerical  simulation. 
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teach.  At  the  downstream  end  of  the  long  pool,  CRM  73.7,  the  wave  shape  and  damping  is  well 
represented  in  the  model.  Only  the  smaller  waves  with  short  wavelengths  are  lagged  in  the  model. 
The  linear  analysis  of  the  model  projected  this  behavior  for  short  period  waves  through  reaches 
where  diffusion  is  large  and  the  value  of  0  required  for  adequate  diffusion  in  the  model  is  small  or 
negative. 

Comparison  between  the  stage  data  at  CRM  71 .4  and  the  numerical  simulation  at  CRM  71 .2  also 
reveals  excellent  model-prototype  agreement.  The  origin  of  the  timing  error  of  the  wave  arriving 
at  hour  136  is  most  likely  cumulative,  as  this  error  is  also  evident  in  the  model-prototype  comparison 
at  CRM  73.7.  Finally,  the  comparison  between  prototype  stage  at  CRM  67.3  and  modeled  stage 
at  CRM  67.2  is  again  excellent.  Due  to  negative  8  values  in  the  model,  small,  short  period  waves 
arriving  between  hours  70  and  90  lag  slightly  in  the  simulation. 

The  success  of  the  simplified  flow  routing  model  of  the  Norris  tailwater  has  important  implica¬ 
tions.  The  presence  of  backwater  reaches  in  natural  rivers  does  not  significantly  affect  unsteady 
flow  waves  and  need  not  be  considered  explicitly.  Pooled  reaches  are  adequately  characterized  by 
a  small  bottom  slope  and  roughness.  Again  in  this  study,  dynamic  waves  were  not  observed  moving 
ahead  of  the  main  flow  or  propagating  upstream  as  a  result  of  a  wave  reflection.  The  reproduction 
of  all  features  of  the  measured  stage-time  traces  by  the  model  demonstrates  that  steep-sided  hydro¬ 
graphs  in  a  mildly  sloping  stream  are  not  affected  significantly  by  inertia.  As  the  relative  importance 
of  the  inertia  terms  is  greatest  for  these  conditions,  it  follows  that  inertia  is  unimportant  in  shallow, 
free-flowing  rivers.  The  extension  of  this  conclusion  to  deeper  rivers  was  indicated  by  Stoker  (1957) 
in  his  study  of  rapidly  rising  flood  waves  on  the  Ohio  River.  He  reported  that  the  first  measurable 
disturbance  traveled  far  behind  the  initial  dynamic  wave  at  the  wave  speed  used  in  simplified 
routing  methods. 

CONCLUSIONS 

Physically  based  relationships  for  the  diffusion  and  dispersion  coefficients  that  describe  wave 
movement  in  rivers  were  developed  from  a  linear  analysis  of  the  dynamic  equations  of  open  chan¬ 
nel  flow.  These  relationships  revealed  that  the  magnitudes  of  the  coefficients  are  related  and  that 
each  is  positive.  By  recasting  the  resulting  equation  in  dimensionless  form,  insight  regarding  the 
relative  magnitudes  of  the  convection,  diffusion  and  dispersion  terms  describing  a  river  flow  wave 
can  be  obtained.  The  magnitudes  of  these  terms  reflect  upon  the  importance  of  the  physical  proc¬ 
esses  affecting  the  flow,  and  provide  justification  for  using  simplified  models.  For  example,  both 
the  adequacy  of  the  kinematic  wave  approximation  and  the  potenital  for  shock  formation  in  the 
channel  are  evaluated  from  the  relative  magnitude  of  the  dimensionless  diffusion  coefficient.  In 
general,  diffusion  of  short-period  waves  is  an  important  process  that  cannot  be  neglected. 

The  diffusion  wave  model  for  flow  in  tailwater  streams  is  a  simplified  inertia-free  routing  model 
that  allows  variable  wave  diffusion  and  does  not  require  a  downstream  boundary  condition.  The 
modified  equation  and  von  Neumann  analyses  provided  stability,  diffusive  and  dispersive  charac¬ 
teristics  of  the  model.  The  analyses  are  complementary,  each  having  its  particular  strengths.  Iden¬ 
tical  stability  conditions,  developed  with  each  approach,  revealed  that  stability  does  not  impose  a 
time  step  limitation  upon  the  model. 

The  dissipative  and  dispersive  characteristics  of  the  model  are  sensitive  to  the  selection  of 
spatial  grid  resolution,  the  Courant  number  and  the  8  parameter.  As  revealed  in  the  von  Neumann 
analysis,  model  damping  increases  as  wavelength  relative  to  the  spatial  grid  length,  Ax,  decreases. 
Both  analyses  showed  that  damping  increases  as  6  decreases  and  as  the  Courant  number  increases. 

An  expression  quantifying  the  numerical  diffusion  of  the  model  was  an  important  result  of  the 
modified  equation  analysis.  This  expression  was  found  to  be  an  accurate  representation  of  nu¬ 
merical  damping.  The  von  Neumann  and  modified  equation  analyses  of  model  phase  error  concurred 


and  were  supported  over  wide  ranges  of  0  and  the  Courant  number  by  the  linear  routing  studies. 

The  frequently  utilized  Hirt  analysis  yielded  incorrect  phase  relationships.  The  routings  and  analyses 
revealed  that  damping  and  phase  errors  occur  when  the  spatial  grid  is  overly  coarse  relative  to  wave¬ 
lengths  of  interest. 

Guided  by  the  linear  analyses,  selection  of  the  spatial  mesh  dimension  and  the  maximum  Courant 
number  can  be  made  a  priori  to  achieve  optima]  accuracy  for  a  given  application.  The  only  calibra¬ 
tion  required  for  model  application  is  the  adjustment  of  Manning’s  roughness. 

The  capabilities  and  accuracy  of  the  diffusion  wave  model  are  enhanced  by  allowing  the  param¬ 
eter  d  to  vary  so  that  a  balance  is  maintained  between  physical  and  numerical  diffusion.  The  numer¬ 
ical  stability  of  the  model  does  not  place  a  lower  bound  on  the  value  of  0.  The  limitation  specified 
in  many  simplified  models-that  the  weighting  parameter  applied  to  the  time  derivative,  in  this  case 
0,  be  greater  than  O.O-should  not  be  generally  applied.  Mildly  sloping  rivers,  modeled  with  a  fine 
spatial  mesh  for  adequate  wave  resolution,  require  negative  values  of  the  weighting  parameter  for 
proper  wave  diffusion. 

A  physical/numerical  dispersion  balance  eliminating  mode)  phase  error  cannot  be  maintained 
in  the  diffusion  wave  model  simultaneously  with  the  diffusion  balance.  Negative  values  of  0  cause 
some  amount  of  lagging  phase  error  in  the  simulation  that  will  be  largest  for  the  shorter  wavelengths. 
The  phase  lag  introduced  in  the  Norris  tailwater  simulation  as  a  result  of  allowing  negative  values 
of  6  was  minimal,  noticeable  only  for  small,  short-period  waves.  Phase  lag  was  not  apparent  in  the 
Apalachia  tailwater  simulation  where  balanced  diffusion  was  maintained  with  positive  values  of  0 . 
For  many  cases  of  interest,  phase  error  can  be  minimized  by  judicious  selection  of  numerical  param¬ 
eters. 

The  release  hydrographs  and  measured  stage  data  from  the  Apalachia  and  Norris  Dam  tailwaters 
provided  a  discriminating  test  of  model  performance  for  wide  ranges  of  bed  channel  slope  and 
roughness.  The  effect  upon  the  flow  of  numerous  large-scale  roughness  elements  in  the  Apalachia 
tailwater  channel  was  adequately  described  by  Manning’s  equation.  The  water  surface  slope  term 
in  the  stage-discharge  relationship  was  necessary  even  in  this  steeply  sloping  river  to  accurately 
model  the  propagation  speed  of  the  flow  releases.  Neglecting  this  term  caused  a  time  lag  of  each 
wave  arrival.  The  inclusion  of  local  inflow  in  the  Apalachia  model  was  needed  to  correctly  propa¬ 
gate  the  small  flow  release.  Correct  wave  celerity  of  this  release  could  not  be  achieved  through 
model  calibration. 

The  ability  of  the  model  to  simulate  flow  in  both  tailwaters  demonstrates  its  generality  for 
rou ring  rapidly  varying  unsteady  flow  in  rivers.  As  expected,  the  release  hydrograph  in  the  Apa- 
lacuij  tailwater  experienced  a  small  amount  of  diffusion  relative  to  that  in  the  Norris  tailwater. 

In  each  case,  the  effect  of  diffusion  upon  the  flow  was  reproduced  with  the  diffusion  wave  model. 
Both  tailwater  studies  indicated  that  backwater  effects  upon  the  flow  of  pooled  river  reaches  need 
not  be  modeled  explicitly. 

The  relative  importance  of  the  inertia  terms  in  the  momentum  equation  is  greatest  for  sleep- 
sided  hydrographs  in  mildly  sloped  streams.  In  the  Norris  Dam  tailwater  test,  however,  dynamic- 
waves  of  a  measurable  size  propagating  in  either  the  upstream  or  downstream  direction  were  not 
observed.  These  observations  c.-upled  with  the  success  of  the  inertia-free  Norris  tailwater  model 
support  analytical  findings  that  inertia  is  not  significant  in  shallow,  free-flowing  rivers  where  flow 
typically  occurs  at  relatively  small  Froude  numbers. 
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